function [xe,ye,e,c] = ellipse(a,b,xc,yc,theta,N)
%
%[xe,ye,e,c] = ellipse(a,b,xc,yc,theta,N)
%
% Draws an ellipse in which
%
%   a      : major half-asis (if undefined ;
%   b      : minor half-axis;
%   (xc,yc): center (if both are unassigned, the default choice xc = yc = 0
%            is assumed);
%   theta  : angle betweem the major axis and the x axis (if undefined, 
%            the default value theta = 0 is considered);
%   N      : number of ellipse points, that is the length of output
%            vectors xe,ye (if undefined, the default N = 100 is used).
%
% The third output is the ellipse eccentricity, and the fourth the focal
% half-distance.
% If a < b, the roles of major and minor axis are inverted.

% G. Teza, 2005, 2021.

if nargin < 6
    N = 100; 
end
if nargin < 5
    theta = 0; 
end
if nargin == 3 
    help ellipse; 
    error('xc and yc must be both defined or both undefined'); 
end
if nargin == 2
    xc = 0; 
    yc = 0; 
end
if nargin < 2
    help ellipse; 
    error('not enough input data'); 
end
if a <= 0 
    help ellipse; 
    error('The major axis must be a positive real number');
end

if a < b
    baux = b; 
    aaux = a; 
    a = baux; 
    b = aaux;
    ind_inv = 1;
else
    ind_inv = 0;
end
c = sqrt(a^2-b^2);
e = c/a;
rhomin = a-c;

phi = linspace(0,2*pi,N);
rho = rhomin*(1+e)./(1+e*cos(phi-theta)); % ellipse in a focus-centered frame 
xer = rho.*cos(phi)+c*cos(theta);
yer = rho.*sin(phi)+c*sin(theta);
if ind_inv == 0
    xerr = xer; 
    yerr = yer;
else
    yerr = xer; 
    xerr = -yer;
end
xe = xerr+xc; 
ye = yerr+yc;